function [foF2, M3000, NmF2] = LegendreF2(mu, lat, lon, CF2, Cm3)
%
% F2-Layer Critical Frequency and M(3000) by Legendre Calculation
%
%DESCRIPTION:
%This function computes the F2-Layer Critical Frequency and M(3000) by
%Legendre Calculation.
%
%PROTOTYPE:
% [foF2, M3000F2] = LegendreF2(mu, lat, lon, CF2, Cm3)
%
%--------------------------------------------------------------------------
% INPUTS:
%   mu         [1x1]       MODIP                     [deg]
%   lat        [1x1]       Latitude                  [deg]
%   lon        [1x1]       Longitude                 [deg]
%   CF2        [76x1]      Coeff.s for foF2          [-]
%   Cm3        [49x1]      Coeff.s for M(3000)F2     [-]
%--------------------------------------------------------------------------
% OUTPUTS:
%   foF2       [1x1]       F2-Layer Crit. Freq.      [MHz]
%   M3000      [1x1]       M(3000)F2                 [-]
%   NmF2       [1x1]       F2-Layer Max. Density     [1e11 m-3]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Constants
DR = pi/180;                %Conversion Factor: deg->rad

%Vectors of Sine and Cosine of Coordinates Memory Allocation
M       = zeros(12,1);
P       = zeros(9,1);
S       = zeros(9,1);
C       = zeros(9,1);
fn      = zeros(9,1);
K       = zeros(9,1);
M3000F2 = zeros(7,1);
H       = zeros(7,1);

%MODIP Coefficients
M(1) = 1;
for k = 2 : 12
    M(k) = sin(mu*DR)^(k-1);
end

%Latitude and Longitude Coefficients
for n = 2 : 9
    P(n) = cos(lat*DR)^(n-1);
    S(n) = sin((n-1)*lon*DR);
    C(n) = cos((n-1)*lon*DR);
end

%F2-Layer Critical Frequency
%Order 0 Term
f1 = 0;
for k = 1 : 12
    f1 = f1 + CF2(k)*M(k);
end
fn(1) = f1;
Q = [12, 12, 9, 5, 2, 1, 1, 1, 1]';
K(1) = -Q(1);
for n = 2 : 9
    K(n) = K(n-1) + 2*Q(n-1);
end
%Higher Order Terms
for n = 2 : 9
    aux = 0;
    for k = 1 : Q(n)
        aux = aux + (CF2(K(n)+2*k-1)*C(n) + CF2(K(n)+2*k)*S(n)) * (M(k)*P(n));
    end
    fn(n) = aux;
end
%Final Summation
foF2 = 0;
for n = 1 : 9
    foF2 = foF2 + fn(n);
end

%Computation of M(3000)F2
%Order 0 Term
M3000F2_1 = 0;
for k = 1 : 7
    M3000F2_1 = M3000F2_1 + Cm3(k)*M(k);
end
M3000F2(1) = M3000F2_1;
R = [7, 8, 6, 3, 2, 1, 1]';
H(1) = -R(1);
for n = 2 : 7
    H(n) = H(n-1) + 2*R(n-1);
end
%Higher Order terms
for n = 2 : 7
    aux = 0;
    for k = 1 : R(n)
        aux = aux + (Cm3(H(n)+2*k-1)*C(n) + Cm3(H(n)+2*k)*S(n)) * (M(k)*P(n));
    end
    M3000F2(n) = aux;
end
%Final Summation
M3000 = 0;
for n = 1 : 7
    M3000 = M3000 + M3000F2(n);
end

%F2-Layer Maximum Density [1e11 m-3]
NmF2 = 0.124 * foF2^2;

end